
########################################

Table A5. 12-Indicators Per Factor

########################################


##############


    N=300

#############

    
library(lavaan)
library(boot)     
    
    
identity<-diag(1, 24,24)
phi<-matrix(c( 1, 0.3, 0.4,
               0.3, 1, 0.5,
               0.4, 0.5, 1
), 3,3)

ld<-matrix(c(0.65, 0.65, 0.7, 0.7, 0.7, 0.7,0.6, 0.5,0.5, 0.5, 0.6, 0.55,
             0.5, 0, 0, 0,  0, 0, 0, 0, 0, 0, 0, 0, 
             0,0,0,0,0,0,0,0,0,0,0,0, 
             
             0,0,0,0,0,0.5,0,0,0,0,0,0,
             0.7, 0.5, 0.5, 0.65,0.5, 0.5,0.6, 0.55, 0.6, 0.45,0.5,0.45,
             0,0,0,0.5,0,0,0,0,0,0,0,0,
             
             0,0,0,0,0,0,0,0,0,0,0,0,
             0,0,0,0.45,0,0,0,0,0,0,0,0,
             0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.7, 0.70, 0.450, 0.5,0.65,0.55
             
), nrow=36, ncol=3)

ld2<-t(ld)
psi<-diag(1,36)-diag(ld %*% phi %*% ld2)
diag(psi)


popmodel<-'
f1 =~ 0.65*x1 + 0.65*x2 + 0.7*x3 + 0.7*x4 + 0.7*x5 + 0.7*x6 + 0.6*x7 
    + 0.5*x8 + 0.5*x9 + 0.5*x10 + 0.6*x11 + 0.55*x12 + 0.5*x13
f2 =~ 0.7*x13 + 0.5*x14 + 0.5*x15 + 0.65*x16 + 0.5*x17 + 0.5*x18
    + 0.6*x19 + 0.55*x20 + 0.6*x21 + 0.45*x22 + 0.5*x23 + 0.45*x24 + 0.5*x6 + 0.5*x28
f3 =~ 0.5*x25 + 0.5*x26 + 0.5*x27 + 0.6*x28 + 0.6*x29 + 0.6*x30 + 0.7*x31
    + 0.7*x32 + 0.45*x33 + 0.5*x34 + 0.65*x35 + 0.55*x36 + 0.45*x16

f3 ~~ 0.4*f1 + 0.5*f2
f1 ~~ 0.3*f2

f1 ~ 0*1
f2 ~ 0*1
f3 ~ 0*1

x1 ~~ 0.5775*x1
x2 ~~ 0.5775*x2
x3 ~~ 0.51*x3
x4 ~~ 0.51*x4
x5 ~~ 0.51*x5

x6 ~~ 0.05*x6
x7 ~~ 0.64*x7
x8 ~~ 0.75*x8
x9 ~~ 0.75*x9
x10 ~~0.75*x10

x11 ~~  0.64*x11
x12 ~~ 0.6975*x12
x13 ~~ 0.055*x13
x14 ~~ 0.75*x14
x15 ~~ 0.75*x15

x16 ~~ 0.0825*x16
x17 ~~ 0.75*x17
x18 ~~ 0.75*x18
x19 ~~ 0.64*x19
x20 ~~ 0.6975*x20

x21 ~~ 0.64*x21
x22 ~~ 0.7975*x22
x23 ~~ 0.75*x23
x24 ~~ 0.7975*x24
x25 ~~ 0.75*x25

x26~~ 0.75*x26
x27~~ 0.75*x27
x28~~ 0.09*x28
x29~~ 0.64*x29
x30~~ 0.64*x30

x31~~ 0.51*x31
x32~~ 0.51*x32
x33~~ 0.7975*x33
x34~~ 0.75*x34
x35~~ 0.5775*x35
x36~~ 0.6975*x36
'

mis_model<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 +x24
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36
f3 ~~ f1
f1 ~~ f2


'



set.seed(312)

mydata<-simulateData(popmodel, model.type = "sem", sample.nobs=300L)
fit1 <- sem(model = mis_model, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")



modindices(fit1, minimum.value = 0, sort = TRUE)[1:12,]


newpar = '
f2	=~	x6
f1	=~	x13
f2	=~	x28
f3	=~	x16
f3	=~	x6
f1	=~	x16
x13	~~	x16
f2	=~	x5
f2	=~	x32
f3	=~	x13
x29	~~	x30
f2	=~	x4'


lavTestScore(fit1, add=newpar)




##############. Bootstrap LM Test


newpar = '
f2=~x6
f1=~x13
f2=~x28
f3=~x16
f3=~x6
f1=~x16
x13~~x16
f2=~x5
f2=~x32
f3=~x13
x29~~x30
f2=~x4'




set.seed(312)
Data<-mydata
# Initialize the matrix to store bootstrap results
bin.1000 <- matrix(NA, nrow = 500, ncol = 2)

# Loop through bootstrap iterations
for (i in 1:500) {
  
  # Perform one bootstrap iteration
  boot.res <- tryCatch({
    # Generate bootstrap sample
    boot.idx <- sample.int(nrow(Data), replace = TRUE)
    Data.boot <- Data[boot.idx, ]
    
    # Fit the SEM model
    fit2 <- sem(model = mis_model, data = Data.boot, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
    
    # Calculate scores
    score1 <- lavTestScore(fit2, add = newpar)$uni[1, 4]
    score2 <- lavTestScore(fit2, add = newpar)$uni[1, 6]
    
    # Return scores
    list(score1 = score1, score2 = score2)
  }, error = function(e) {
    # If SEM model did not converge, return NULL
    NULL
  })
  
  # Check if boot.res is not NULL (i.e., SEM model converged)
  if (!is.null(boot.res)) {
    # Update the matrix with scores
    bin.1000[i, 1] <- boot.res$score1
    bin.1000[i, 2] <- boot.res$score2
  }
}

# Calculate the means of the scores
t1 <- mean(bin.1000[, 1], na.rm = TRUE)
t2 <- mean(bin.1000[, 2], na.rm = TRUE)
t1
t2





#######  Bootstrap W test

w_model<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + b2*x13 + b6*x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 +x24 + b1*x6 + b3*x28 + b8*x5 + b9*x32 + b12*x4
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + b4*x16 + b5*x6 + b10*x13
f3 ~~ f1
f1 ~~ f2

x13~~b7*x16
x29~~b11*x30

'



fit2 <- sem(model = w_model, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


con='

b1==0'



set.seed(312)
Data<-mydata
# Initialize the matrix to store bootstrap results
bin.1000 <- matrix(NA, nrow = 500, ncol = 2)

# Loop through bootstrap iterations
for (i in 1:500) {
  
  # Perform one bootstrap iteration
  boot.res <- tryCatch({
    # Generate bootstrap sample
    boot.idx <- sample.int(nrow(Data), replace = TRUE)
    Data.boot <- Data[boot.idx, ]
    
    # Fit the SEM model
    fit2 <- sem(model = w_model, data = Data.boot, meanstructure = FALSE, likelihood = "wishart", estimator = "ML")
    
    # Calculate Wald test statistics and p-values
    wald_stat <- lavTestWald(fit2, constraints = con)$stat[1]
    wald_p_value <- lavTestWald(fit2, constraints = con)$p.value[1]
    
    # Return statistics
    list(stat = wald_stat, p_value = wald_p_value)
  }, error = function(e) {
    # If SEM model did not converge, return NULL
    NULL
  })
  
  # Check if boot.res is not NULL (i.e., SEM model converged)
  if (!is.null(boot.res)) {
    # Update the matrix with Wald test statistics and p-values
    bin.1000[i, 1] <- boot.res$stat
    bin.1000[i, 2] <- boot.res$p_value
  }
}

# Calculate the means of the Wald test statistics and p-values
w1 <- mean(bin.1000[, 1], na.rm = TRUE)
w2 <- mean(bin.1000[, 2], na.rm = TRUE)
w1
w2



###############################################################

      Likelihood Ratio Test

###############################################################

m_0<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36
f3 ~~ f1
f1 ~~ f2

'
lrt_0<-sem(model = m_0, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")




m_1<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36
f3 ~~ f1
f1 ~~ f2

'
lrt_1<-sem(model = m_1, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


m_2<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36
f3 ~~ f1
f1 ~~ f2

'
lrt_2<-sem(model = m_2, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")



m_3<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36
f3 ~~ f1
f1 ~~ f2

'
lrt_3<-sem(model = m_3, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


m_4<-'  


f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16
f3 ~~ f1
f1 ~~ f2


'
lrt_4<-sem(model = m_4, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


m_5<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6
f3 ~~ f1
f1 ~~ f2
'
lrt_5<-sem(model = m_5, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")



m_6<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6
f3 ~~ f1
f1 ~~ f2

'
lrt_6<-sem(model = m_6, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")



m_7<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x6 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6
f3 ~~ f1
f1 ~~ f2

x13~~x16

'
lrt_7<-sem(model = m_7, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


m_8<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x6 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28 + x5
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6
f3 ~~ f1
f1 ~~ f2

x13~~x16
'


lrt_8<-sem(model = m_8, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


m_9<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x6 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28 + x5 + x32
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6
f3 ~~ f1
f1 ~~ f2
x13~~x16
'


lrt_9<-sem(model = m_9, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")



m_10<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x6 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28 + x5 + x32
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6 + x13
f3 ~~ f1
f1 ~~ f2
x13~~x16
'

lrt_10<-sem(model = m_10, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


m_11<-' 
f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x6 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28 + x5 + x32
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6 + x13
f3 ~~ f1
f1 ~~ f2

x13~~x16
x29~~x30
 '
lrt_11<-sem(model = m_11, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")



m_12<-' 
f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x6 + x16
f2 =~ x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x6 + x28 + x5 + x32 + x4
f3 =~ x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x16 + x6 + x13
f3 ~~ f1
f1 ~~ f2

x13~~x16
x29~~x30
 '

lrt_12<-sem(model = m_12, data=mydata, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")


lavTestLRT(lrt_0, lrt_1, lrt_2, lrt_3, lrt_4, lrt_5, lrt_6, lrt_7, lrt_8, lrt_9, lrt_10, lrt_11, lrt_12)





########.  The End

####### All results match the manuscript, 12/4/24
